Determination of the chemical potential using energy-biased sampling 
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An energy-biased method to evaluate ensemble averages requiring test-particle insertion is pre- 
sented. The method is based on biasing the sampling within the subdomains of the test-particle 
configurational space with energies smaller than a given value freely assigned. These energy-wells are 
located via unbiased random insertion over the whole configurational space and are sampled using 
the so called Hit&Run algorithm, which uniformly samples compact regions of any shape immersed 
in a space of arbitrary dimensions. Because the bias is defined in terms of the energy landscape it 
can be exactly corrected to obtain the unbiased distribution. The test-particle energy distribution is 
then combined with the Bennett relation for the evaluation of the chemical potential. We apply this 
protocol to a system with relatively small probability of low-energy test-particle insertion, liquid 
argon at high density and low temperature, and show that the energy-biased Bennett method is 
around five times more efficient than the standard Bennett method. A similar performance gain is 
observed in the reconstruction of the energy distribution. 
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I. INTRODUCTION 

The chemical potential is a central quantity under- 
pinning many physical and chemical processes, such 
as phase equilibria, osmosis, thermodynamic stability, 
binging affinity and so on». However, its evaluation 
by computer simulation is more complicated and time- 
consuming than for other intensive thermodynamic quan- 
tities, such as the pressure P or temperature T. While 
P and T can be evaluated from averages over mechanical 
properties of molecules (forces, velocities and positions), 
the chemical potential is a thermal average and therefore 
it requires sampling the phase space of the system. In- 
deed, computing the chemical potential is a special case 
of the more general problem of computing a free-energy 
difference A\ — Aq between two states (labelled as and 
1), a problem for which the inherent difficulty is well 
understood- j2 i 3 i 4 . Free energy perturbation (FEP) is an 
important category of methods for free energy calcula- 
tion; we refer to the recent works by Lu et alX and by 
Shirts and Pande^ for review and comparisons. As ex- 
plained by Lu et alJL, the general working equation for 
FEP methods can be cast as 



cxp[-/3(A 1 - A )} = 



(w(u) exp[-/3u/2]) 
(ifl(u)expH3ti/2]}i' 



(1) 



with (3 = 1/ksT and u = U\ — Uq the energy differ- 
ence between both systems; Kb is the Boltzmann con- 
stant. The angular brackets denote ensemble averages 
performed on the system labelled by the subscript "0" or 
"1" . The weighting function w(u) is arbitrary and differs 
for each method introduced in the literature. 

The chemical potential is the free energy difference be- 
tween two thermodynamic states differing by the pres- 
ence of a single molecule. In other words, the chemical 
potential is A\ - A where Ai = A(N + l,V,T) and 
A = A(N, V, T). Here A(N, V, T) is the Helmholtz free 



energy of the system which depends on the number of 
molecules N, the volume V and temperature of the sys- 
tem. In order to express the averages of Eq. Q in terms 
of one-dimensional integrals of the energy difference u one 
can then introduce the following distribution functions 6 



f(u) = / {5{u-U l + U )) (i V- 1 dv, 



g(u) = (6{u-Ui + U ))i 



(2) 
(3) 



where <5(.) is the Dirac delta function. In Eq. 
U\ = ?7i(R , r), where is the configuration of the 
first N molecules and r denotes the configuration of the 
N+l molecule. Note that in Eq. (0) the N+l molecule 
acts as a "test-molecule" which probes the system "0" 
(i.e. the system with N molecules), but does no inter- 
act with it. Therefore f(u) is the probability density of 
the N molecule ensemble increasing in potential energy 
by an amount u if this test-molecule were randomly in- 
serted into the ensemble. Conversely, g(u) is the proba- 
bility density of the (N+l)-molecule ensemble decreasing 
in potential energy by an amount u if a randomly selected 
real molecule were removed from the ensemble. 

^From Eq. an expression for the excess chem- 

ical potential /i = A\ — Aq — fiid (where \iid is the ideal 
gas chemical potential^) can be derived in terms of the / 
and g distributionsi£*£ 



exp(Pn) 



J w(u)g(u)du 



J w(u)f(u) exp(— j3u)du 



(4) 



A good choice of the weighting function w[u) is key for 
the efficiency of the method. For instance, the Widom 
method^ (w(u) = 1) is known to provide very poor con- 
vergence at large densities. The Widom method is a sin- 
gle stage FEP, meaning that sampling is only performed 
in the reference system "0" (i.e., in the / distribution, 
see Eq. (@J). As discussed by Lu et aim , multiple staging 
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provides much better efficiency. The efficiency is gener- 
ally defined as the reciprocal of the product of the vari- 
ance of the estimator multiplied by its cost n cos t (that is, 
the total number of energy evaluations performed by the 
algorithm) 



£ = (n cost Var[/?/i] 



(5) 



Bennett^ showed that the variance of Eq. (0J is min- 
imised if the weighting function is w(u) — T[f3{u — c)], 
where T{x) = l/(f + exp(a;)) is the Fermi function and c 
is an arbitrary constant. The Bennett estimator is then 



/3fi = In 



f (T[-P(u-c)]) g 
\ {T[(3{u - e)]) f 



(3c, 



(6) 



where the subscripts g and / indicate (simple) averages 
over the distributions g(u) and f(u). The value of c pro- 
viding the minimum variance and maximum overlap is 
c = [i and to evaluate fj, using the optimum c(= fi) 
one requires to use a self-consistent procedure, iterat- 
ing the value of c in Eq. © and resetting c = fi un- 
til {T[—(3{u — c)]) g = {T[fl(u — c)])f. In practise, this 
step only requires a small number of iterations. Recent 
publications 1,5 demonstrate that the Bennett method re- 
mains the best general method to compute the chemical 
potential for many applications. 

Note that the Bennett method is a two-stage FEP and 
therefore it also requires sampling of the system "1" . In 
the case of the determination of the chemical potential 
this system has N+l molecules and g(u) is obtained from 
its single-molecule energy distribution. However this ex- 
tra requirement is not really a drawback. Lu et alX 
showed that, provided N > 0(100), the g— average can 
be evaluated in the same simulation as is used to sample 
the / distribution (system "0") without any noticeable 
loss in accuracy. The g distribution (constructed from 
the energy of the real particles) is thus a byproduct of 
the simulation so the average {T) g does not demand any 
extra computational cost. 

Another group of methods for determination of the 
chemical potential are based on biased instead of uni- 
form sampling. In particular, cavity-biased methods first 
select spherical cavities of minimum radius R c (a free pa- 
rameter) in which to insert the test-molecule. This ac- 
celerates the evaluation of the ensemble average in dense 
phases because the low-energy configurations of the test- 
molecule (with large Boltzmann factors) are usually lo- 
cated in larger cavities with less steric hindrance. Vari- 
ations of this method have been proposed by several au- 
thors; these include the Cavity Insertion Widom method 
(CIW) due to Mezei and coworkers**, the Excluded Vol- 
ume Map Sampling by Deitrick et al& and the method 
proposed by Pohorille and WilsoniS. The cavities are 
located by a grid search over the whole simulation cell. 
A cavity centre is assigned at each grid point whose dis- 
tance to the closest particle is greater than R c . In order 
to correct the bias introduced in sampling only inside the 
cavities one also has to calculate the probability of find- 
ing a cavity, which is obtained in the same grid-search 



step. A drawback of the cavity-biased method is that it 
is only indirectly related to the test-particle energy via 
the excluded volume. This fact introduces a certain inac- 
curacy in the estimation of the chemical potential, as it 
can depend on the value of the cavity radius R c selected. 
For instance, the CIW has recently been used to calcu- 
late the chemical potential of several species across a lipid 
bilayer«. As a test calculation the authors estimated the 
chemical potential of water in water and reported varia- 
tions of about 1 Kcal/mol as R c was varied from 2. 6 A to 
2.8A Also, using R c £ [2.6, 2.9]A resulted in uncertain- 
ties of about 2 Kcal/mol in estimates of the excess chem- 
ical potential of some species across the lipid layer. Note 
that the important region of the cavity-biased method is 
constructed over the translational degrees of freedom of 
a "coarse-grained" spherical molecule with an effective 
radius. This means that it can only be applied to small 
solutes with spherical or roughly spherical shapes^. 

In this work we present an energy-biased method for 
the estimation of the chemical potential and reconstruc- 
tion of the energy distribution f(u) in dense phases. The 
idea is to restrict the sample to an important region 
defined by the set of bounded domains in the configu- 
rational space of the test-molecule where the energy u 
is smaller than a given free parameter u w . We denote 
as an energy-well each compact subdomain within the 
test- molecule energy-landscape for which u < u w . Note 
that the present approach retains the main benefit of 
the cavity-biased method, but provides an exact evalu- 
ation of the energy distribution f(u) and the chemical 
potential, because the energy-wells arc defined directly 
in terms of the energy landscape. Moreover our energy- 
biased method does not assume any particular molecu- 
lar shape and therefore it may be used for non-spherical 
molecules and can coherently sample over rotational de- 
grees of freedom as well. 

We also note that the number of stages are not lim- 
ited to two. When systems and 1 are very different it 
may be impossible within the simulation time to sample 
the importance region of the two systems. In this case 
it is more efficient to compute the total free energy dif- 
ference by using a set of intermediate states. The energy 
bias method can be applied on each of these intermediate 
state transitions at the cost of performing independent 
simulations for each state. Other approaches include, for 
instance, slow and fast growth methods where the sys- 
tem is changed from one state to another within a cer- 
tain simulation time r (large for slow growth). The fast 
growth method consists of sampling rapid transformation 
from many simulations which are then combined by us- 
ing Jarzynski nonequilibrium work relationii to obtain 
the total free energy difference. 

The rest of the paper proceeds as follows. The energy- 
biased method is explained in Sec. |n] while in Sec. IIIII 
we derive an analytical expression for the efficiency of the 
method and estimate the optimal parameter u w by max- 
imising the efficiency. In Sec IIVI the method is tested 
in liquid argon at high density (modelled as Lennard- 
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Jones atoms) where it is used to reconstruct the test- 
particle energy distribution f(u) and the chemical po- 
tential. We also demonstrate the gain in efficiency ob- 
tained with energy-biased sampling with respect to uni- 
form sampling. We conclude with a summary of our 
findings in Sec. E] Finally in Appendix E] we briefly 
explain the Hit&Run algorithm which efficiently samples 
bounded regions of arbitrary shape immersed in an arbi- 
trary number of dimensions. 



II. OVERVIEW OF THE METHOD 

As stated in the introduction, energy-biased sampling 
consists of uniform sampling of the importance region 
defined by the set of subdomains in the test-molecule 
configurational space where its potential energy is less 
than u w . The probability density is therefore given by 



h(u) 



f(u)/F w u<u n 







u > u v 



(7) 



where the normalisation factor F w = f™™ f{u)du is the 
cumulative probability of the unbiased distribution f(u) 
and u w is an arbitrary energy (free parameter). 

Note that the energy-biased distribution of Eq. (7J can 
be straightforwardly combined with any of the popular 
methods to calculate the chemical potential from Eq. 10}. 
We shall use the Bennett method due to its excellent 
performance. Introducing the weighting function w(u) = 
!F\P(c — u)] in Eq. JBJ and using Eq. (0, one obtains the 
energy-biased Bennett estimator for 



Pfj, = In 



F w {F c ) f 



(3c, 



(8) 



where we have introduced the notation T c = T[(3{u — cj] 
to indicate that after the ensemble average we still have 
a function of c. As before, the subscript h indicates the 
average over the biased distribution of Eq. (Q) . 

Sampling from the energy probability distribution h(u) 
requires a more careful consideration of the energy land- 
scape of the system. We indicate by r a configuration of 
the (N+l)th molecule and by R the configuration of the 
remaining N molecules. For a simple argon fluid r£fl 
where D C R 3 , while for a 3 sites flexible water model like 
TIP3P D C i? 9 , which includes the three Euler angles de- 
termining the molecule orientation, the H-O-H angle and 
the two H-0 distances. 

As shown in Fig. the region 



A*. 



{r e D : u(r, R) < u w } 



(9) 



is composed of many disconnected bounded regions of 
different sizes such that A Uw — il a A"^, where each 
is now a connected region. Of course, for u w — > oo we 
have that all the regions A" connect and A^ = D, the 
entire domain. The sampling algorithm must reproduce 




FIG. 1: Energy landscape for the three-dimensional configu- 
rational space generated by inserting an argon atom in a cube 
of side 16 A of argon fluid. The isosurfaces of regions A UqlI 
are shown for u w equal to 1 (dark grey) and to 10 Kcal/mol 
(light grey). 



a uniform probability distribution 



Pu m (r) 



fl(A Uti 



(10) 



where VL(A Um ) is the volume of the region. 

For a given energy bias u w , the algorithm for selecting 
configurations r according to Eq. i|l(J|) can be described 
in terms of two main steps which are applied iteratively: 

1. Locate a compact energy- well A" in the configu- 
rational space D, where u < u w . 



2. Sample the energy- well A" 
ability density. 



with a uniform prob- 



The simplest procedure for locating energy wells in 
step (1) is to perform a random search over the whole 
configurational space until a fixed number of cavities is 
found. This procedure, however, does not avoid the prob- 
ability of exploring the same well more than once, and 
we observed that it can easily lead to highly correlated 
data. Instead we perform step (1) by choosing points 
on a grid within the whole configurational space of the 
test-molecule. In the case of the Lennard- Jones fluid, the 
three-dimensional configurational space is probed at the 
nodes of a Cartesian grid of size n x x n y x n z , where 
n a is the number of nodes along the coordinate a. We 
observed that the minimum distance between nodes that 
guarantees statistically independent samples is around 
0.5ct. 
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An energy well is found at each node where the energy 
of the test- molecule is u < u w . Then, the locations of 
each of these nodes are used as starting configurations 
for independent well samplings. In this way we ensure 
that we are sampling different cavities for each explored 
configuration (snapshot) of the system. Note that using 
grid-sampling the number of cavities found per snapshot 
is a fluctuating quantity. 

The search requires an average of n = 1/F W energy 
evaluations to locate one well (i.e. one configuration with 
energy u < u w .) During this same step (I) one can cal- 
culate the cumulative probability F w from the estima- 
tor m/riQ, with no being the total number of samples 
(Bernoulli trials) and m the number of successful trials 
with u < u w , i.e., the total number of energy- wells found. 
This number m/no converges to F w as no — > oo and, for 
a finite number of statistically independent trials no, its 
variance is (I — F w )F w /riQ. In practise, the estimation 
of F w requires the number of unbiased samples to be 
no >> \/F w ; this condition also ensures that a signifi- 
cant number of energy-wells (m > 0) are to be found. 

Step (2) of the loop mentioned above requires a pro- 
cedure to sample in an unbiased way the interior of each 
energy well. This is a delicate step because any bias in- 
curred in sampling the importance region will be trans- 
fered to the estimator for resulting in inaccuracy of 
the method. To tackle this problem we use the so-called 
Hit&Run algorithm 12 , which is explained in Appendix 

El 



III. EFFICIENCY AND OPTIMAL 
PARAMETERS OF THE METHOD 

We now calculate the efficiency of the method and pro- 
vide a way of choosing the optimal value of the parameter 
u w by maximising the efficiency. We also compare the 
efficiency of the estimator in Eq. © based on energy- 
biased sampling with that of the standard Bennett algo- 
rithm of Eq. ©. 



A. Energy-biased Bennett method 



Let us now consider the variance of the estimator in Eq. 
(JSJ) , which is the sum of the variance of the estimator for 
F w and the estimator for the ensemble average 



Var EB [/3 M ] =Var[lnF ro ]- 



1 



1 



1 



u(^F c )h noF w n w {T c )h 
(13) 

where we have used the relation Var[ln(F lu )] ~ 
Va.r[F w ]/Fl = (1 - F w )/(n F w ) ~ l/(n F w ), for F w « 
1. Here no is the number of random insertions in the 
entire configurational space and n w is the number of in- 
dependent samples within the importance region u < u w . 

The probability of finding an energy-well with u < u w 
using uniform sampling over the whole configurational 
space is F w , so the number of cavities found after no trials 
is ttl — F w hq. If the number of statistically independent 
samples per well is s, the total number of independent 
samples within the restricted configurational space u < 



u w is 



n w = n sF w 



(14) 



We note that the number of independent samples per well 
s depends on the fluid considered and, of course, on the 
biasing energy u w . In Appendix iBl we provide a way of 
estimating s from the outcome of the data obtained from 
Hit&Run sampling. Inserting Eq. (|14|l into Eq. Ijl3(l one 
obtains for the energy-biased algorithm 



Va.r EB [13 n] 



1 

no 



1 



(15) 



In deriving Ea. fT5|l we used that (T c ) / = F w {T c )h up 
to a negligible amount. This can be seen by noticing 
that the function JF[/3(u — c)] in the integrand of (J- c ) / = 
J^oo f( u )3~[(3(u — c)]du decays exponentially for u > c. 
Hence, in any practical case (u w > c) most of the integral 
weight comes from u < u w , for which the energy-biased 
reconstruction of the energy profile f(u) is exact (see Fig. 

El- 

We now evaluate the cost, which is given by the total 
number of energy evaluations of the test molecule needed 
to obtain n w samples: 



The variance of the Bennett method can be cast in 
terms of the probability densities f(u) and g(u). Starting 
from Eq. ©, after some algebra the variance of the 
Bennett method assumes the form 



Var B [/3/x] = 



1 



noW(«-c)])/' 



(11) 



where no is the number of insertions used to sample 
the complete configurational space of the test-particle. 
Note that the computational cost of the standard Ben- 
nett method is no, so according to @ and Eq. i|llfl its 
maximum efficiency is given by 



o/f- 



(12) 



n cos t = n + n^/a, 



(16) 



where a < 1 is the acceptance ratio of the Hit&Run 
sampling algorithm, defined in Appendix ^ Introduc- 
ing Eq. (fT^f) into Eq. ljPo|l we obtain 



n-cost = n 1 + 



sF, 



(17) 



For the energy-biased algorithm the efficiency is e = 
(ncostVarss^u])" 1 . Using Eq.JTU) and Eq.(H3) one ob- 
tains 



£ eb — 



1 

J- in 



1 s 

s{F c ) f a 



F 

■*■ 11. 



a (^"c)j 



(18) 
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By maximising the efficiency e = e(F w ) in Eq. I|18|) with 
respect to F w , one obtains the optimal value F^ pt and 
the maximum efficiency e_EB max = £EB(F° pt ): 



popt _ 



' EB„ 



s{T c )f' 



(19) 
(20) 



Finally we compare the efficiency of the energy-biased 
algorithm with that provided by the Bennett algorithm, 
given by Eb = (J~ c )f. According to Eq. (|20|l the ratio of 
efficiencies is given by 



SB 



£ EB„ 



= 2 



{?c)f . s(f c )f 



(21) 



Equation (|21|l yields the range of values of (F c ) j 
for which the energy-biased Bennett estimator for 0fi 
method is more efficient than the standard (unbiased) 
Bennett algorithm. Note that for s = -J a/ {T c ) / the effi- 
ciency ratio given by Eq. I|21|l reaches its minimum value, 
£fl/£BB m „ = 4-^/ (Tc) //a, and therefore sb < £eb if 
{T c ) f > a/16. Hence the energy-biased method is suited 
for fluids at high densities or low temperatures or for 
molecular fluids with low insertion probability. In this 
regime (J- C )f << a/16 and the dominant term in Eq. 
(|21|l is 1/s, hence e_EB max — ssb- In other words, the 
maximal efficiency of the present energy-biased method 
is limited by the average number s of independent sam- 
ples that can be obtained within one energy-well. As 
shown in Appendix [5] for the Lennard-Jones fluid we 
have observed that in the most unfavourable case (high 
density and low temperature) s ~ [5 — 10]. 



Using Var(if) = H(l - H)/n w and Var(F w ) = F W (F W 
1)/ no one obtains 



Var BB 



F W (F W — \)H 2 
n 



H(l - H)Fl , F W H(1 
1 



H) 



n w n Q n w 

(23) 

Note that, as expected, for H ~ 1 one recovers the vari- 
ance of the unbiased insertion method. The interesting 
part of the energy distribution is the importance region, 
located in the low energy range, where H << 1. In this 
regime one can make the approximation 1 — H ~ 1 . Using 
F = HF W and n w — uqsF w , one gets 



1 



F F 1 
Var BS = —— + --, 

n \F W s n sF u 



(24) 



Note that the term in brackets is the reduction in vari- 
ance with respect to uniform unbiased sampling. Because 
F w is evaluated from no probes, this means that necessar- 
ily n Q » l/F w so the third term inside the brackets is 
much smaller than unity. On the other hand, for the low 
energy range considered F << F w and one finally con- 
cludes that Var^s ~ Var(F)/s, where Var(F) ~ F/n is 
the variance obtained in the unbiased uniform sampling 
of the whole domain. 

The cost associated with the energy-biased procedure 
is n cost = n (l + sF w /a). In the case of a Lennard- 
Jones liquid we have found that a ~ 0.17 and s ~ O(10), 
while the optimal cumulative probability is F w < 10~ 3 . 
This means that, in practical situations, sF w /a < 1 and 
n cos t > no- Thus, according to Eq. (22} the energy- 
biased sampling procedure is around s times faster than 
a uniform unbiased (grid or random) sampler in recon- 
structing the low energy range of f(u). As before, s is the 
average number of independent samples taken per well. 



B. Reconstruction of the energy distribution 



We now show that the reconstruction of /(it) using the 
energy-biased procedure (EB) is faster and more efficient 
than that obtained using any unbiased sampler which 
uniformly explores the whole configurational space. To 
that end we consider the evaluation of the cumulative 
probability F(u) = /_ f(u')du for u < u w (i.e. for 
F(u) < F w ). We shall compare the variance of two es- 
timators for F: one based on uniform insertion over the 
whole domain and the other based on the energy-biased 
procedure. The variance of the unbiased estimator is sim- 
ply Var(F) = F(l — F)/n and for low energies (F « 1) 
its efficiency is l/F. The expected value of the energy- 
biased estimator is HF W , where H(u) = f_ h(u')du' 
is the cumulative probability of the biased distribution 
in Eq. Q. This estimator is constructed as a product 
of two statistically independent fluctuating variables and 
its variance isi^ 

Var BB (F) = Var(HF w ) = F*Vax(H) (22) 
+ i? 2 Var(^) +Var(F w )Var(if). 



IV. RESULTS 

In order to confirm the foregoing theoretical relations 
about efficiency and variance reduction, we performed 
molecular dynamics simulations of a Lennard-Jones liq- 
uid at high density and low temperature (p — 0.0236A" 3 
and T — 84K). These simulations were performed in a 
cubic periodic box of side L = 10a. We used the stan- 
dard Verlet method^ to integrate Newton's equations of 
motion, incorporating a Langevin thermostat^ to keep 
the system in the NVT ensemble. 

During the simulation, the iterative loop (l)+(2) ex- 
plained in Sec. [HI was performed m times per time in- 
terval Stsamp = 0.5t, which corresponds to about three 
times the collision time. The search for wells performed 
in step (1) was done by probing at the nodes of a Carte- 
sian grid comprising 15 3 nodes. This ensured that the 
explored cavities are independent. All the cavities found 
in step (1) were sampled using the Hit&Run algorithm 
(see Appendix lAl . 
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A. Estimation of the chemical potential 

One way to measure the efficiency of the method 
is to evaluate the convergence of the estimated value 
of the chemical potential for an increasing number of 
test-particle probes n cost . Convergence can be calcu- 
lated from the difference between successive values of p n , 
where n(= n cost ) indicates the total number of evalua- 
tions of the test-particle energy. Figure [21 shows how 
this difference decreases in calculations based on both 
the energy-biased and the unbiased samples. These cal- 
culations correspond to liquid argon with number density 
p = 0.0236A -3 and temperature T = 84K (these values 
correspond to p — 0.92cr 3 and T = 0.7 in Lennard-Jones 
units), for which the average of the Fermi function is 
{T c ) f = 8.9 x 10~ 6 . According to Eq. JTJJJ the optimum 
value of F w is 0.0012, which corresponds to u w ~ 14.19 
Kcal/mol. We selected the predicted optimum parameter 
(u w = 14.19 Kcal/mol) and performed d = 15 samples 
per well. As can be seen in Fig. |2 for equal numbers 
of energy probes ( tl — ii'cost): the average difference be- 
tween successive estimates of the chemical potential via 
the energy-biased method is about five times smaller than 
that obtained with the unbiased sampler. As predicted 
by Eq. J2U, such a gain in efficiency is consistent with 
the average number s of independent samples per well 
(see table Ej , which for this simulation was s ~ 5. 

Evaluations of the chemical potential for Lennard- 
Jones (LJ) fluids are shown in Table [I] together with 
the estimated efficiency of each calculation. For a LJ 
fluid with p = 0.02360l~ 3 and T = 84K the numeri- 
cally obtained net gain is around 7, which coincides with 
the prediction in Eq. I|21|) using s — 7. For illustra- 
tive purposes we also analysed a case for which the effi- 
ciency of our implementation of the energy-biased sam- 
pling is similar to the uniform-unbiased Bennett method. 
For instance, (F c ) f = 0.0102 for p = 0.01755A~ 3 and 
T = 178. 5K. Using a = 0.165 and the (optimum) num- 
ber of samples s = */ a/ (JF C ) j ~ 4 in Eq. (j21(l one ob- 
tains Eb / £eb^ x — 1; our numerical calculations, with 
u w = 7.33 and d — 8, confirmed this conclusion. We note 
that for any value of u w considered the energy-biased es- 
timation of the chemical potential p agrees within about 
0.01 Kcal/mol with the unbiased Bennett result. This is 
illustrated in Table ITT1 where we show the estimated p for 
the higher density liquid, using several values of u w . 
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FIG. 2: (Top) The estimation of the chemical potential 
plotted against the overall number of energy probes (n cost ). 
We compare the standard (uniform sampling) Bennet and 
Widom methods with the corresponding energy-biased ver- 
sions of these methods. The calculations correspond to a 
Lennard-Jones fluid with p = 0.0236A" 3 , and T = 84K 
(p = 0.92<t 3 and T = 0.7e/k B in reduced LJ units); the 
energy-biased sampling was done using u w = 14.19 Kcal/mol 
and d = 15 samples per well. (Bottom) The convergence 
measured as the squared difference between consecutive es- 
timations with increasing cost (An = 10 s ). In the energy- 
biased method the cost is given by n C03t = no (1 +dF w /a), 
where no is the number of random probes used to evaluate 
F w (— F(u w ) — 0.00122) and a = 0.165 is the acceptance ra- 
tio of the Hit&Run sampler. Circles are the successive values 
of (p n — Pn+An) 2 (shown only for the energy-biased case) and 
solid lines are the average over twenty consecutive differences. 



B. Reconstruction of the energy distribution /(it) 

In Fig. |2| we compare the reconstructed energy dis- 
tribution f(u) at energies u < u w with that computed 
from an unbiased method, which consists of a large num- 
ber of random insertions within the entire configurational 
space. Figure 13 clearly illustrates that the energy-biased 
method exactly reproduces the unbiased distribution f(u) 
for energies smaller that u w . This attractive feature is a 
consequence of the fact that it is easy to exactly cor- 



rect for the bias in terms of the cavity energies. This is 
not true for the accessible volume of the molecule, as in 
cavity-biased procedures^. 

In order to illustrate the above conclusion we show 
in Fig. 21 the estimation of the cumulative probabil- 
ity F(u) versus the total number of test-particle en- 
ergy probes used for the evaluation. The particular case 
shown corresponds to u — 5 Kcal/mol, for a LJ liquid at 
p = 0.02361" 3 and T = 84K. The energy-biased sam- 
pling was done using u w = 14.19 Kcal/mol and d = 15 
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T(K) 


[iEB 


/is 






seb/eb 


F w 


d 


0.02360 


84 


-0.336 


-0.323 


0.9 x 10 5 


(1.2 ±0.1) x 10 4 


7± 1 


0.00122 


15 



TABLE I: Comparison of the chemical potential (in Kcal/mol) calculated via the standard Bennett method (i.e. using uniform 
unbiased sampling) and the energy-biased Bennett (subscript EB). The inefficiency of both methods (reciprocal of efficiency) is 
also shown. In the case of the standard Bennett method we write the minimum inefficiency (e^ 1 = (•?>)/) while the inefficiency 
of the energy-biased method was obtained from numerical calculation of the variance of /3fi, using block-analysis (see Appendix 
B or e.g. Refsp-*-) and agrees within error bars with the theoretical expression of Eq. 1181 (see text). The error in eeb/sb 
comes mainly from the uncertainty in the numerical calculation of Var^ b ■ 
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FIG. 3: The energy distribution f(u) obtained from 4 x 10 6 
random insertions over the whole configurational domain is 
compared with energy-biased sampling in the restricted con- 
figurational space u < u w . The calculations correspond to 
the same case as in Fig. |5] The energy cavities are sampled 
using the Hit&Run algorithm, which provides an unbiased 
reconstruction of the energy distribution for any value of u w 
chosen. 



FIG. 4: The cumulative probability F(u) = /"^ f{u')du' for 
u = 5 Kcal/mol versus the total number of energy evaluations 
of the test-particle n coa t . The liquid is the same as in Fig. [5] 
We compare the estimations of F(u) for grid-sampling (with 
a regular mesh of 36 3 nodes) and for energy-biased samplings 
within u < u w — 14.19 Kcal/mol, performing d — 15 samples 
per well. The cumulative probability at u w is F w = F(u w ) = 
0.00122. 



samples per well, and for this calculation we obtained 
s ~ 5 (see Appendix and Table HTjl. Compared with 
the unbiased procedure, the reduction of variance pro- 
vided by the energy-biased sampler is immediately ap- 
parent on inspection of Fig. 0] A numerical evalua- 
tion of the variance of each data set in Fig. 0] pro- 
vides: Var^s = 4.14 x 10~ 5 /rio, while the (best) result 
for the algorithm based on uniform unbiased sampling is 
Va.r(F) = F/no — 1.9 x lCT 4 /n . Hence the net gain 
in efficiency is about 4.6, in agreement with the value of 
s = 5 obtained from the independent correlation analy- 
sis explained in Appendix [51 As shown in Table [I] the 
estimated net gain in the evaluation of the chemical po- 
tential compared with the unbiased Bennett method is 
7 ± 1, which is close to the estimate s ~ 5 obtained from 
the analysis of the cumulative probability. 



V. CONCLUSION 

We have presented a new method for sampling the 
energy of a test-molecule in order to calculate single- 
particle ensemble averages and, in particular, the chem- 
ical potential. The method, called energy-biased sam- 
pling, restricts the important region to the bounded do- 
mains in the test-molecule energy-landscape where the 
test-molecule energy u is smaller than a given free param- 
eter u w . This energy-biased sampling retains the princi- 
pal benefit of cavity-biased methods^ in the sense that, 
by sampling only within regions with a significant Boltz- 
mann factor, convergence is greatly accelerated with re- 
spect to uniform sampling. Furthermore, because the 
energy-biased sampling is accurately defined in terms of 
the test-particle energy it has some important benefits: 
first, it allows accurate reproduction of the test-particle 
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energy distribution f(u) and the chemical potential; sec- 
ond, it is possible to sample cavities of arbitrary shape 
(not only spherical ones) and to generalise the cavity 
dimensionality to include the rotational degrees of free- 
dom in the energy-well reconstruction; finally, and rather 
importantly, it enables one to combine the sampling re- 
sults with standard free energy perturbation (FEP) for- 
mulae. In particular, we combined it with the Bennett 
method 8 which minimises the variance of the estimator 
and has proved to be the best method in the literature^. 
Energy-biased sampling is a general protocol to bias the 
sampling and consists of two sequential steps: (1) search- 
ing and (2) sampling the interior of energy-wells. In this 
work we have implemented these two steps using rela- 
tively simple algorithms: uniform unbiased search and 
Hit&Run sampling. However we note that other solu- 
tions are also possible. For instance, non-uniform sam- 
pling of the importance region may surely increase the 
efficiency of the present method. In dense systems, the 
searching step becomes the most difficult one and a more 
effective extension of this method could be to perform 
a biased search (using, for instance, some variation of 
the USHER algorithmic*^) so as to significantly increase 
the probability of finding favourable cavities for insertion 
of the test particle. These extensions are left for future 
studies. 
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i. Choose a random direction e and find the intersec- 
tions of the cavity border with the line r(A) = ro+Ae, 
where A is a real number. As the cavity A is bounded 
the intersection is composed by two points r(A + ) and 
r(A") (here A+ > and A~ < 0). 

ii. Select a point ri within the segment (r(A + ), r(A - )), 
i.e., 

ri =r(A-)+£(r(A+)-r(A-)) (A2) 

where £ 6 (0, 1) is a uniformly distributed random 
number. 

iii. Sample at rx, set ri — > ro as the new starting point 
and go to (i). 

The above procedure is repeated to obtain the desired 
number of samples d. In our case the starting point 
for the sample chain, r , is the test-particle configura- 
tion returned by the algorithm for energy-well searching 
(£/(r ,R) < u w ). In order to locate the borders of the 
energy well r(A + ) and r(A - ) we use the following proce- 
dure. Starting from ro we cross the well along the line 
defined by the random unit vector e moving in steps of 
size 5s, i.e., according to 

r(fc) = r + kSse, (A3) 

with k being an integer starting from k = ±1. The en- 
ergy is computed at each point r(k) until one crosses 
the edges of the well at k = k + and k = k~ (for which 
w(r(k ± ), R) > u w ). An approximate location of the cav- 
ity borders is provided by setting X^ 1 = k^. We used 
typically Ss ~ 0.3A and required, on average, about five 
iterations to cross the well in one random direction (this 
value depends on the density and u w ). Note that the 
acceptance ratio is a = (fc + — fc - ) -1 and for the high 
density cases considered here a ~ 0.17. 



APPENDIX A: SAMPLING BOUNDED REGIONS 
WITH THE HIT&RUN ALGORITHM 



APPENDIX B: OPTIMAL NUMBER OF 
SAMPLING DIRECTIONS 



There exists a relatively large literature on sampling 
a bounded connected region (sec for instance Ref.— and 
references therein). In this work we have used the so- 
called Hit&Run algorithm for its simplicity and good 
performance^. The Hit&Run sampler is a special Monte 
Carlo Markov chain which draws numbers from an as- 
signed distributioni2*ii p(r), where r G A lies within 
a bounded connected region of an n-dimensional space 
A C R n . In our case, p(r) is a uniform probability den- 
sity over the region A"^ such that 



p(r) 



1 



(Al) 



The Hit&Run algorithm starts from a point ro within 
the bounded region A and performs the following steps: 



It is possible to reduce the cost without increasing the 
variance by setting the number of samples per cavity d 
equal to or somewhat larger than s, the average number 
of independent samples per cavity. Note that the number 
of statistically independent samples within one cavity is 
s = d/r c , where r c is an empirically estimated autocor- 
relation length of the whole chain of data. This number 
t c can be estimated from the large m limit of the quan- 
tity 77iVar[.F( m) ]/Var[.F], where T c = F\p{u - c)] is the 
Fermi function evaluated at a single energy u and T^™^ 
denotes the mean of m consecutive T values. 

The value of s can be estimated by performing sev- 
eral Hit&Run samplings with an increasing number of 
directions per cavity d > s, then computing r c for the 
chain of samples and evaluating d/r c , which should be 
nearly independent of d. We carried out this evaluation 
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of s for varying values of u w within the same system 
and for fixed u w and varying density. The results of this 
study, reported in Table EI clearly indicate that s does 
not greatly vary for a broad range of values of the cavity- 
border energy u w . In fact, at low and moderate values 
of u w the energy-cavities are isolated and their average 
size (in A) grows quite slowly with u w . This is due to 
the steepness of the hard-core part of the Lennard-Jones 
potential. Above a certain energy u w the cavities be- 
come connected and a steep rise in the average size of 
the energy-cavities is observed. This is reflected in the 



value of s. As shown in Table UTI for u w = 14.19 Kcal/mol 
we obtained s ~ 4.5 and s ~ 11 for two calculations us- 
ing d = 15 and d = 100 respectively. We obtained a 
relatively close value s ~ 7 for twice as large an energy 
limit u w = 28.38 Kcal/mol. However using u w = 165.53 
Kcal/mol the average number of independent samples in- 
creased up to 25, reflecting the more complex shape and 
larger volume of these energy cavities. In summary, for 
the optimum range of values of u w ~ [10 — 30] Kcal/mol 
we find s ~ [5 — 10] in the case of the Lennard-Jones 
liquid. 
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u w (Kcal/mol) 


d 


no 




s 


ricost/no 


H (Kcal/mol) 


28.38 


20 


4.2 10 5 


1.1 10" 3 


7 


1.5 


-0.32 


28.38 


100 


2.95 10 6 


1.08 10" 3 


7 


3.8 


-0.353 


14.19 


100 


4.3 10 6 


2.66 10~ 3 


12 


1.7 


-0.335 


14.19 


15 


1.0 10 7 


2.77 10~ 3 


5 


1.1 


-0.334 


165.53 


200 


2.76 10 5 


1.3 10" 5 


25 


85.8 


-0.357 



TABLE II: Details of the energy-biased calculations in a Lennard- Jones (LJ) liquid at density p = 0.0236A -3 and temperature 
T = 84K (p — 0.92 and T — 0.7 in LJ units). We compare the results for varying values of the energy parameter u m , samples 
per cavity d and varying number no of energy probes within the unbiased distribution. The cumulative probabilities of the 
unbiased distribution (F w = F(u w ) = J^f(u)du) are F„(165.53) = 0.0704, F w (28.38) = 0.00458, F w ( 14. 19) = 0.00122. The 
average of the Fermi function in the biased distribution (J- C )h is defined in Eq. JSJl. The average number of independent 
samples per cavity s is obtained from a = d/r c , where the correlation number r c is calculated from the correlation between the 
whole chain of data. The overall number of energy probes in the energy-biased method is n cost = no (1 + dF w /a), where a is 
the acceptance ratio obtained for the Hit&Run sampler, a ~ 0.17. The estimation of the chemical potential using the standard 
(unbiased) Bennett method with 1.1 x 10 7 energy samples is u = —0.323 Kcal/mol. 



